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Abstract — A  new  parametric  inversion  technique  to  estimate  a 
buried  pipe  location  was  developed  for  borehole  radar  cross-hole 
measurements.  This  technique  evaluates  the  shapes  of  the 
approximated  and  measured  arrival  time  curves  instead  of  the 
first  arrival  time  itself  or  wavefield  in  conventional  inversion  or 
tomographic  techniques.  In  this  study,  we  propose  an  algorithm 
of  the  technique  and  demonstrate  its  performance  for  synthetic 
and  measured  datasets. 
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I.  Introduction 

To  determine  locations  of  buried  pipes  such  as  gas  and 
water  supplies,  and  subsurface  cavities  is  important,  for 
example  before  starting  excavations.  For  borehole  radar,  the 
horizontal  resolution  is  relatively  poor  and  it  is  difficult  to 
determine  the  horizontal  location  of  pipes  with  high  accuracy. 
Tomography  and  inversion  techniques  have  generally  been 
used  for  such  investigations  [1],  [2].  Moreover,  migration 
based  imaging  techniques  also  have  been  proposed  [3]. 
However,  these  methods  require  densely  sampled  datasets  and 
are  computationally  expensive.  It  is  often  difficult  to  acquire 
appropriate  datasets  for  those  techniques  due  to  temporal  and 
spatial  limitations  of  the  surveys  especially  in  urban  areas.  In 
this  contribution,  we  propose  a  new  inversion  scheme  for 
bistatic  transmission  radar  measurements. 

Conventional  radar  tomography  techniques  try  to  image 
contrasts  of  electrical  medium  parameters.  In  other  inversion 
techniques,  e.g.,  [1]  and  [2],  the  total  wavefield  is  modeled  and 
compared  to  the  measured  wavefield.  In  contrast,  the  new 
inversion  technique  to  be  described  here  compares  the  shape  of 
the  first  arrival  time  curves  of  measured  data  with  that  of 
theoretical  curves.  Furthermore,  our  inversion  technique  can 
use  known  information  efficiently  and  explicitly.  Unlike  other 
inversion  techniques,  the  result  of  this  technique  is  not  an 
image,  but  an  explicit  indication  of  the  pipe  location. 

This  paper  describes  the  algorithm  of  the  new  inversion 
technique.  Then  it  is  applied  to  synthetic  datasets  and  measured 
datasets  to  estimate  a  buried  metallic  pipe  location.  For  the 
synthetic  datasets  by  finite  difference  time  domain  (FDTD) 
method,  the  technique  is  demonstrated  for  a  homogeneous  and 
heterogeneous  media  to  show  the  robustness  of  this  technique 
against  the  inhomogeneity. 
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II.  Algorithm  of  the  Inversion  Technique 

A.  Calculation  of  Approximated  Arrival  Times 

Suppose  that  the  target  whose  location  we  are  trying  to 
estimate  is  a  metallic  pipe  with  a  known  diameter  buried  in  a 
homogeneous  medium.  Two  vertical  boreholes  are  located  in 
the  plane  perpendicular  to  the  pipe.  Then  approximately,  the 
first  arrivals  are  given  by  signals  that  propagate  in  this  plane 
where  the  propagation  paths  are  the  shortest.  The  shapes  of  the 
first  arrival  time  curves  resulting  from  those  paths  are  different 
for  various  locations  of  the  pipe,  transmitter,  and  receiver.  If 
the  pipe  is  located  between  the  antennas,  the  signal  cannot  pass 
through  the  metallic  pipe  and  propagates  along  the  curved 
surface  of  the  pipe  as  shown  in  Fig.  1,  thus  delaying  the  first 
arrival  time.  In  this  case  the  approximated  first  arrival  time 
tcai(z,  x,  zh  zr)  can  be  calculated  using  the  ray-path  model  with  a 
propagation  velocity  in  the  medium  v,  an  assumed  pipe 
location  (z,  x),  and  the  approximated  path  length  L  for  the 
transmitter  at  zt  and  the  receiver  at  zr  as  follows. 

teal  V,  x,z„zr)  =  L(z,  x,zt,zr)/v  ( 1 ) 

B.  Evaluation  of  Curve  Shapes 

The  approximated  first  arrival  curves  are  compared  with  the 
arrival  time  curves  from  the  measured  datasets.  Similarities  of 
those  curves  are  calculated  for  various  assumed  pipe  locations, 
and  a  location  that  shows  the  highest  similarity  indicates  the 
true  location  of  the  pipe.  Picking  first  arrival  times  is  difficult 


Figure  1 .  Model  for  calculation  of  the  approximated  first  arrival  time. 
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because  of  noise  and  limited  frequency  bandwidth.  Generally, 
inversion  techniques  that  use  the  first  arrival  times  strongly 
depend  on  the  quality  of  the  pickings.  Here  we  propose  two 
schemes  for  evaluating  the  first  arrival  curve  shapes  that  do  not 
depend  on  the  manual  pickings. 

1)  Gradient  Error  Scheme :  It  is  easy  to  pick  the  times 
tmeas(z t,  zr)  at  the  maximum  or  a  thresholded  amplitude  of  each 
trace.  The  curves  obtained  in  this  way  are  no  true  first  arrival 
time  curves,  but  we  assume  them  to  be  parallel  to  the  true  first 
arrival  time  curves.  Then  similarity  in  shape  between  the 
picked  curve  tmeas  and  the  modeled  curves  tcai  can  be  evaluated 
by  taking  the  error  e(z,  x,  zt)  of  the  gradients  of  tmeas  and  tcah 
which  is  given  by  integration  along  the  receiver  depth. 

e(z, x,z,)  =  j dzr  |3z  tmeas (z, ,zr)~ dz  tcal (z, x, z, , zr )|  (2) 

This  error  can  be  calculated  for  each  transmitter  position.  The 
total  error  for  a  series  of  transmitter  positions  is  given  by 
physical  conjunction  of  the  errors  e(z,  x,  zt)  and  can  be 
calculated  by  integration  along  the  transmitter  depth. 

eto,al  (z> x)  =  j  dz,  e{z,x,zt)  (3) 

In  this  scheme,  the  lower  error  values  indicate  a  higher 
similarity.  Therefore,  the  minimum  error  value  indicates  a 
location  that  has  the  highest  possibility  of  the  true  target 
location. 

2)  Correlation  Scheme :  From  the  modeled  first  arrival 
time  curve,  we  firstly  construct  a  binary  reference  H(z ,  x,  zt9  zr, 
t).  The  reference  is  a  matrix  equal  in  size  to  the  measured 
datasets.  The  elements  of  H  that  correspond  to  the 
approximated  first  arrival  times  are  set  to  1  and  the  others  are 
set  to  0.  Cross-correlation  between  this  reference  and  the 
measured  dataset  S(zt,  zn  t)  indicates  a  similarity  of  the  first 
arrival  time  curves,  and  the  similarity  is  defined  by  the 
following  expression. 

R(z,x,zt )  =  max|j*J dzdzr  S(zt,zr,T)H(z,x,zt,zr,t -r)|  (4) 

To  obtain  the  total  correlation,  such  similarity  values  are  to  be 
integrated  along  the  transmitter  depth. 

K,al(Z’X)  =  \dz,R(Z’X’Z,)  (5) 

In  this  scheme,  the  higher  correlation  values  indicate  a  higher 


similarity.  Therefore,  the  maximum  correlation  value  indicates 
a  location  that  has  the  highest  possibility  of  the  true  target 
location.  Comparing  gradient  error  scheme,  the  similarities  can 
be  evaluated  without  any  picking  operations  by  this  scheme. 

By  using  the  two  described  schemes,  we  can  avoid  the  error 
that  occurs  for  picking  the  first  arrival  times  from  a  measured 
dataset,  and  we  do  not  need  to  consider  how  to  define  them. 
This  makes  our  inversion  scheme  simple,  robust  and  easy  to 
operate. 

III.  Application  to  Synthetic  Datasets 

The  inversion  technique  was  applied  to  synthetic  datasets 
generated  by  3-D  FDTD  simulation  of  cross-hole  fan 
measurements  for  the  configuration  shown  in  Fig.  2.  The 
distance  between  two  boreholes  is  3.3  m,  and  a  1.0  m  diameter 
metallic  pipe  is  buried  at  a  depth  of  12.0  m  and  at  a  distance  of 
2.0  m  from  the  transmitter  borehole.  The  depth  of  the 
transmitter  ranged  from  ll.Omto  13. 0  m  with  a  0.5  m  step  and 
of  the  receiver  from  10.0  m  to  14.0  m  with  a  0.1  m  step. 

A.  Homogeneous  Media  Model 

For  a  homogeneous  media  model,  we  set  a  permittivity  of 
20 Cq  and  a  conductivity  of  0.01  S/m  to  the  whole  surrounding 
media.  Fig.  3  shows  the  gradient  error  maps  for  each 
transmitter  depth.  The  total  error  map  obtained  by  integration 
of  all  the  error  maps  is  shown  in  Fig.  4(a),  and  the  total 
correlation  obtained  from  all  the  correlation  maps  is  shown  in 
Fig.  4(b).  In  the  error  and  correlation  maps  for  a  single 
transmitter  position,  the  low  error  and  the  high  correlation 
regions  are  not  localized  (see  Fig.  3).  In  contrast,  the  total  error 
and  correlation  maps  have  a  localized  region  as  shown  in  Figs. 
4(a)  and  (b).  The  result  from  both  evaluation  schemes  indicates 
a  location  that  is  exactly  the  same  as  the  modeled  location  of 
the  metallic  pipe. 


Figure  2.  FDTD  Model. 
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Figure  3.  Gradient  error  maps  for  the  trasmitter  at  depths  of  (a)  1 1.0  m,  (b)  1 1.5  m,  (c)  12.0  m,  (d)  12.5  m,  and  (e)  13.0  m. 
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Figure  4.  The  results  of  the  inversion  technique  calculated  from  synthetic 
datasets  with  homogeneous  media  model  using  (a)  the  gradient  error  scheme 
and  (b)  the  correlation  scheme.  The  crosses  indicate  the  minimun  value  in  (a) 
and  the  maximum  value  in  (b). 
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Figure  6.  The  results  of  the  inversion  technique  calculated  from  synthetic 
datasets  with  heterogeneous  media  model  using  (a)  the  gradient  error  scheme 
and  (b)  the  correlation  scheme.  The  crosses  indicate  the  minimun  value  in  (a) 
and  the  maximum  value  in  (b). 


Position  x  [m]  Relative  Permittivity 

Figure  5.  Relative  permittivity  distribution  of  the  simulated  heterogeneous 
media,  (a)  Vertical  slice  aty  =  5.0  m.  (b)  Histogram  of  the  simulated  relative 
permittivity. 


Figure  7.  Geometrical  sketch  of  the  site. 


B.  Heterogeneous  Media  Model 

For  demonstrating  a  more  realistic  situation,  a 
heterogeneous  model  was  simulated  in  3-D  space  employing  a 
stochastic  fractal  model  [4].  The  simulated  heterogeneous 
media  model  and  its  histogram  are  shown  in  Fig.  5.  This  model 
consists  of  14  permittivities  ranging  from  16.56b  to  23.06b  and  a 
constant  conductivity  of  0.01  S/m.  This  distribution  has  the 
mean  permittivity  of  19.56b  and  the  standard  deviation  of  1.1 6b- 
In  the  calculation  of  the  inversion  technique,  we  set  the  single 
permittivity  of  the  medium  to  20.06b  assuming  homogeneous 
media.  Fig.  6  shows  the  total  error  map  and  the  total  correlation 
map.  Despite  the  localized  regions  are  deformed  comparing 
with  the  homogeneous  model  case  shown  in  Fig.  4,  they  have 
the  minimum  error  and  the  maximum  correlation  at  the  exact 
modeled  location. 

IV.  Application  to  Measured  Datasets 

Field  measurements  were  carried  out  in  an  urban  area  in 
Sendai,  Japan,  to  determine  the  location  of  a  metallic  water 
pipe.  The  pipe  is  0.9  m  in  diameter  and  the  borehole  separation 
is  3.3  m  as  shown  in  Fig.  7.  In  cross-hole  fan  measurements  at 


this  site,  the  transmitter  was  set  in  the  borehole  BH1  at  depths 
of  11.0-13.0  m  with  a  0.5  m  step,  and  the  receiver  was  set  in 
the  borehole  BH2  at  10.0-14.0  m  with  a  0.1  m  step.  An 
approximate  depth  of  the  pipe  and  the  propagation  velocity  in 
the  surrounding  media  were  known  from  previous 
measurements.  Fig.  8  shows  the  raw  radar  profiles  acquired  at 
each  transmitter  depth.  In  these  data,  the  influence  of 
geological  discontinuities  and  inhomogeneity  is  visible, 
especially  in  Figs.  8(d)  and  (e).  To  prevent  it  from  obscuring 
the  inversion  results,  some  of  the  traces  are  excluded.  Wile  the 
error  and  the  correlation  maps  calculated  for  each  transmitter 
position  are  not  localized  as  in  the  case  for  the  synthetic 
datasets,  the  total  error  and  correlation  maps  have  a  localized 
region  as  shown  in  Figs.  9(a)  and  (b).  The  minimum  error  in 
the  total  error  map  is  obtained  at  a  depth  of  12.1  m  and  at  a 
distance  of  2.1  m  from  the  borehole  BH1,  and  the  maximum 
correlation  in  the  total  correlation  map  is  obtained  at  a  depth  of 
12.0  m  and  at  the  distance  of  2.1  m  from  BH1.  These  results 
are  consistent  with  known  information  from  other 
measurements  that  the  pipe  is  located  at  a  depth  of  about  12.0 
m  and  the  borehole  BH2  is  closer  to  the  pipe  than  BH1 . 


Figure  8.  Radar  profiles  obtained  by  cross-hole  fan  measurements  with  the  transmitter  at  (a)  1 1 .0  m,  (b)  1 1 .5  m,  (c)  12.0  m,  (d)  12.5  m,  and  (e)  13.0  m  depth. 
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Figure  9.  The  results  of  the  inversion  technique  for  measured  datasets  using 
(a)  the  gradient  error  scheme  and  (b)  the  correlation  scheme.  The  crosses 
indicate  the  minimum  value  in  (a)  and  the  maximum  value  in  (b). 
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Figure  10.  Velocity  distribution  obtained  by  travel  time  tomography. 


V.  Discussion 

By  using  the  proposed  two  schemes  to  evaluate  arrival  time 
curve  similarities,  we  could  avoid  the  subjectivity  of  picking 
the  true  first  arrival  times.  Moreover,  we  do  not  have  to  specify 
so  many  parameters.  It  needs  to  be  specified  where  we  pick 
times  representing  first  arrival  for  the  gradient  error  scheme, 
while  for  the  correlation  scheme,  what  kind  of  the  reference 
used.  Here  the  reason  why  the  binary  reference,  which  consists 
of  only  0  and  1,  was  used  is  to  reduce  the  parameters  and 
consequently  subjectivity. 

This  inversion  technique  works  well  to  estimate  a  buried 
pipe  location.  However,  for  measured  datasets,  the  estimated 
locations  in  the  both  of  two  total  maps  are  not  localized  as  clear 
as  for  the  synthetic  datasets.  It  is  caused  by  the  transmitter 
positions  where  the  datasets  were  acquired.  The  maps  for  each 
transmitter  depth  are  spreading  along  the  lines  connecting  the 
transmitter  positions  and  the  center  of  the  buried  pipe  as  shown 
in  Fig  3,  and  we  cannot  estimate  the  buried  pipe  location  from 
those  maps.  Whereas  the  total  maps  can  have  localized  region 
by  taking  a  summation  of  all  the  maps,  therefore,  the  datasets 
acquired  by  the  transmitter  at  various  position  especially  both 
above  and  bellow  the  pipe  are  necessary.  In  our  case,  the 
datasets  acquired  by  the  transmitter  bellow  the  buried  pipe 
cannot  be  used  due  to  the  strong  effect  of  the  geology  as  we 
can  observe  in  Figs.  8(d)  and  (e).  That  is  the  reason  why  the 
maps  spread.  However,  it  could  estimate  the  pipe  location, 
although  the  datasets  are  few  and  sparse,  and  conversely,  it 
might  be  possible  to  determine  a  buried  location  using  only  two 
datasets,  which  acquired  at  above  and  bellow  the  target. 

A  travel  time  tomography  for  these  measured  datasets  could 
not  image  the  pipe  clearly  as  shown  in  Fig.  10,  which  is  due  to 
sparse  sampling  and  the  narrow  measurement  region.  Therefore, 
this  inversion  technique  is  more  suitable  to  determine  a  buried 
pipe  location  than  tomographic  technique  for  those  conditions. 

VI.  CONCLUTION 

The  new  inversion  technique  for  estimation  of  a  buried  pipe 
location  is  proposed.  For  the  synthetic  datasets  with 
homogeneous  media,  the  inversion  technique  could  retrieve  the 
exact  location  of  the  buried  metallic  pipe  using  both  the 
gradient  error  and  correlation  scheme,  also  for  the  datasets  with 


the  simulated  heterogeneous  media  model.  This  indicates  the 
inversion  technique  is  rather  robust  against  the  inhomogeneity. 
Moreover,  it  could  estimate  the  pipe  location  for  the  measured 
datasets  in  agreement  with  the  known  information.  This  shows 
that  our  technique  works  better  than  tomography  for  pipe 
location  from  sparse  and  small  datasets.  Further,  it  only 
requires  rather  simple  calculations  that  can  be  done  quickly. 
Therefore,  this  technique  is  robust,  simple,  stable,  and  easy  to 
operate.  It  means  it  can  easily  be  applied  on-site,  which  is  an 
enormous  advantage  for  the  practical  use  of  this  technique  for 
underground  surveying. 

We  think  that  our  inversion  technique  can  be  applied  to 
other  targets,  for  instance  tunnel-like  cavities,  with  small 
modifications.  Because  the  technique  uses  shape  variations  of 
the  arrival  time  curves,  and  it  can  be  applied  to  a  target  that 
delays  or  quickens  the  propagation  time.  Moreover,  it  can  be 
applied  to  other  bistatic  radar  transmission  measurements  such 
as  vertical  radar  profiling  (VRP)  measurements  by  changing 
the  geometry.  Therefore,  it  has  a  high  potential  for  application 
to  various  targets  and  measurements. 
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